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Abstract. We use the Matsubara functional renormalization group (FRG) to 
i describe electronic correlations within the single impurity Anderson model. In contrast 

£h to standard FRG calculations, we account for the frequency-dependence of the two- 

i 1 particle vertex in order to address finite-energy properties (e.g, spectral functions). 

gH ■ By comparing with data obtained from the numerical renormalization group (NRG) 

Q ' framework, the FRG approximation is shown to work well for arbitrary parameters 

, O , ■ (particularly finite temperatures) provided that the electron-electron interaction U is 

not too large. We demonstrate that aspects of (large U) Kondo physics which are 
04 ' described well by a simpler frequency-independent truncation scheme are no longer 

captured by the 'higher-order' frequency-dependent approximation. In contrast, at 
\ small to intermediate U the results obtained by the more elaborate scheme agree 

(Nl ■ better with NRG data. We suggest to parametrize the two-particle vertex not by 

three independent energy variables but by introducing three functions each of a single 



o 



' frequency. This considerably reduces the numerical effort to integrate the FRG flow 

■ equations. 

o 



PACS numbers: 71.27.+a, 73.21.La 
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1. Introduction 

Great advances in nanotechnology over the past few years led to the fabrication and 
intense experimental study of low-dimensional electron systems (quantum dots and 
wires). From a long-term perspective, such systems are of interest as they give insight 
into the physics of more complicated nanodevices needed for quantum information 
processing [1]. The smallness of quantum dots leads to fairly large energy level spacings 
and at sufficiently low temperatures only a few levels are relevant for the description of 
the physics. The latter is then strongly affected by the repulsive interaction between 
the electrons, manifesting e.g. in Coulomb blockade behavior [2] or Kondo screening 
[3]. The physics of quantum wires is described by the Luttinger rather than the Fermi 
liquid theory, in clear contrast to most bulk materials where the effects of Coulomb 
correlations do not extend beyond mild renormalization of Fermi liquid parameters. [4]. 

From the theoretical point of view, a many-particle method is needed to properly 
account for the electron-electron interaction U, which is a vital ingredient to any model 
one might devise to describe the experimental setups. A perturbative approach works 
in some cases, but it fails to describe strong correlations (such as Kondo physics [5]). 
Moreover, it breaks down completely for certain classes of low-dimensional systems 
because of infrared divergences in low-order Feynman diagrams [5] . The latter motivates 
the application of renormalization group (RG) based methods which successively address 
all energy scales of the system, starting from high energies where infrared singularities 
are cut out. One particular implementation of Wilson's general RG idea [6] is the 
numerical renormalization group (NRG), which was originally developed to address the 
Kondo model [7] but later on applied to various kinds of quantum impurity systems 
[8]. It provides a very reliable tool to investigate physical properties of models with 
Coulomb interaction at low energies. However, its applicability is practically limited to 
small systems with a few correlated degrees of freedom because of the computational 
resources required. 

A different RG based approach to solve the quantum many-particle problem is the 
functional renormalization group (FRG) which exactly describes all vertex functions of 
the corresponding system in terms of an infinite hierarchy of coupled flow equations [9]. 
In practice, this hierarchy needs to be truncated (usually by neglecting the flow of the 
three-particle vertex), rendering FRG an approximate method. Moreover, in the context 
of quantum dots and wires the frequency dependence of the flowing two-particle vertex 
was discarded, alltogether resulting in a closed finite set of coupled differential equations 
whose numerical solution gives renormalized frequency- independent system parameters 
embodying the effects of the two-particle interaction [10]. This approximation scheme 
was successfully applied to describe Luttinger liquid behavior of one-dimensional 
quantum wires with local inhomogeneities [12, 13, 14], whereas ordinary perturbation 
theory in the two-particle interaction U is plagued by infrared singularities. In contrast, 
perturbation theory is usually regular for quantum dot models (such as the single 
impurity Anderson model) and there is no inherent need for an RG based framework. 
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However, the appearance of the exponentially (in U) small Kondo energy scale T% 
motivates RG resummations of certain classes of diagrams. The application of the 
frequency-independent functional renormalization group to single- as well as multi-level 
spinful and spin-polarized quantum dot geometries turned out to give surprisingly good 
results in the strong coupling limit, even though the approximation can a priori be 
justified only for small to intermediate U [15, 16, 17, 18]. The most striking observation 
is that the FRG describes aspects of Kondo physics (e.g., an exponential energy scale) 
contained within the single impurity Anderson model (SIAM) very accurately [16]. The 
frequency-independent approximation can thus be viewed as a kind of RG-enhanced 
Hartree-Fock theory which does not suffer from typical mean-field artifacts (such as the 
breaking of spin symmetry). 

By construction, a truncation procedure which disregards all frequency 
dependencies cannot be expected to give reliable results for finite-energy properties. 
Indeed, the SIAM linear-response conductance at zero temperature T = (which is 
a zero-energy quantity) is described well by this level of approximation [16], whereas 
the conductance at T > (which is a finite-energy property) is not. In addition, finite 
energy effects become important if one is interested in extending the method to the non- 
equilibrium situation. It is thus very reasonable to devise a truncation procedure which 
includes the frequency-dependence of the two-particle vertex in order to describe finite- 
energy properties. A first step in this direction was done in Ref. [20] where results for 
the SIAM were presented, illustrating that such a generalization is in principle possible 
[21]. However, there was no systematic study of all system parameters (in particular 
finite temperatures and finite magnetic fields) and the question whether strong-coupling 
physics is captured (as it is partially by the frequency-independent approximation) was 
not answered conclusively. Finally, technical details about how the flow equations are 
actually implemented numerically were not elaborated on. 

It is the first goal of this paper to systematically study the SIAM using a frequency- 
dependent FRG scheme, benchmarking this approximation against data we obtain from 
the numerical RG framework. The second objective is to present technical details 
of this generalization in order to address the issue of numerical artifacts originating 
from the discretization of the Matsubara frequency axis. Even though little is to be 
learned about the physics, we think that systematically discussing the details of this 
very natural and obvious generalization of a frequency-independent FRG truncation 
scheme in its application to the SIAM is of importance, particularly if one aims at 
treating more complex (multi- level or multi-impurity) quantum dot systems. In short, 
the strong-coupling behavior extracted from the FRG approximation turns out to be 
worse than one might have expected, particularly if one has in mind the success of 
the frequency-independent approach. However, at small to intermediate U/T < 5 (r 
being the impurity-lead hybridization) the agreement with NRG reference data improves 
if the more elaborate finite-frequency FRG scheme is employed. The same holds 
if the two-particle vertex is parametrized not by three independent energy variables 
but (numerically far less demanding) by three functions each depending on a single 
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argument. Since the computational effort in solving the flow equations grows only as 
a power law (and not exponentially) with the number of impurities and channels and 
since there is no numerical need to stick to special (symmetric) system parameters, the 
frequency-dependent FRG approximation presented in this paper can be regarded as a 
fast and reliable tool to describe the small- to intermediate-coupling physics of correlated 
quantum dot models which cannot be treated using the NRG. 

The paper is organized as follows. In Sec. 2, the single impurity Anderson model 
is introduced briefly. Next (Sec. 3.1), we recapitulate the general idea of the Matsubara 
functional renormalization group. The flow equations for the SIAM are derived explicitly 
in Sec. 3.2 and Appendix A. We present a cutoff procedure which allows for treating 
both the zero-temperature limit as well as finite T > (Sec. 3.3). Numerical issues 
which come along with the need for discretizing frequency space are commented on 
in Sec. 3.4. Finally, we introduce three functions each of a single energy variable 
to parametrize the two-particle vertex in order to speed up numerics (Sec. 3.5). In 
Sec. 4.1, spectral functions obtained from the various FRG schemes are benchmarked 
against NRG reference calculations. We discuss the average impurity occupation and 
the zero-temperature spectral weight at the chemical potential in Sec. 4.2 and try to 
extract the Kondo temperature from the FRG formalism by considering effective masses 
and static spin susceptibilities (Sec. 4.3). An outlook is given in Sec. 5. 



2. Model 



The Hamiltonian of the single impurity Anderson model [23] consists of three parts: 

H = -ffbath + H dot + if coup- (1) 

The bath is described by noninteracting electrons with a single-particle dispersion 

H hath = J2^k c l a c k^ ( 2 ) 

ha 

— * 

where c^ a denote fermionic annihilation operators for electrons with momentum k and 
spin direction a. The impurity Hamiltonian contains both single- and two-particle terms, 




with d a annihilating an electron of spin a located at the impurity. The single-particle 
energy e was shifted such that e = corresponds to the point of particle-hole symmetry, 
and U and B denote the strength of the Coulomb repulsion and of the magnetic field, 
respectively. Finally, the coupling between the dot and the bath is given by 

#coup = ~tJ2 C t d - + HX - > ( 4 ) 

(T 

where c a is the local electron operator c a = J2^c^ a /\^N. 
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In order to apply the FRG scheme presented in the next Section to the SIAM 
Hamiltonian Eq. (1), the noninteracting bath has to be integrated out using a standard 
projection technique [24]. Thereafter, instead of dealing with an infinite system one only 
needs to consider two interacting (spin up and down) particles. For the noninteracting 
impurity Green function we obtain 

G ° {ZUJ) = ^- e -aB/2 + ,sgnMr ' ^ 

where the hybridization T = ii\t\ 2 J2k^(^ ~ e %)l^ is assumed to be energy-independent 
(wide-band limit). 



3. Method 



In this Section we briefly present the general idea of the Matsubara functional 
renormalization group and discuss the technical details of its application to the SIAM. 
We employ a truncation scheme in which the frequency- dependence of the two-particle 
vertex is kept and introduce an infrared cutoff which allows for treating both T = and 
finite temperatures. The corresponding flow equations are given explicitly in Sec. 3.2 
and Appendix A. We elaborate on the issue of numerical artifacts originating from the 
discretization of frequency space. Further approximations (particularly the frequency- 
independent FRG scheme) are discussed. The numerical renormalization group (which 
will serve to benchmark all FRG results) and its application to the SIAM have been 
presented extensively in the literature (for a review see Ref. [8]), and we will refrain 
from commenting on this method in the present paper. 



3.1. FRG - general idea 

The FRG is one implementation of the general renormalization group idea for interacting 
quantum many-particle systems [9]. It starts with introducing a cutoff A into the 
noninteracting Green function Q°(iu). Here, we choose an infrared cutoff in Matsubara 
frequency space, 

g°(iu) — > G°> A (iuj) = e A (M - A) g°(iu)), (6) 

where the form of Ba is to be specified later on (Sec. 3.3). By virtue of this replacement, 
all irreducible m-particle vertex functions 7 A acquire a A-dependence, and differentiating 
each 7 A with respect to A yields an infinite hierarchy of flow equations which can be 
computed straight-forwardly using functional integrals [9, 25]. In principle, integrating 
from A = 00 (where all 7 m are known) down to A = yields exact expressions for all 
vertex functions. Since this is practically impossible, a truncation procedure has to be 
devised (see below). 

The exact equation for the self energy S A = 7 A reads [9] 



«9 A £ A (l'; 1) = -t£)S a (2; 2')7 2 A (1'2'; 12), 

22' 



(7) 
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with the initial condition being £ A =°° = —JJ /2. The arguments are a shorthand for a set 
of Matsubara frequencies and arbitrary single-particle quantum numbers. The quantity 
<S A denotes the (single-scale) propagator, 



S 1 



(8) 



where Q is the full interacting Green function for a given A, 



(0O,A)-1 _ s a' 



(9) 



Apart from the self-energy itself, the two-particle vertex 7 A couples into Eq. (7), so we 
need to consider its flow equation as well. One can show that it is given by [9] 



9 A72 A (1'2'; 12) = ~TJ2 5 A (3;3')^ A (4;4'){ 72 A (3'4'; 12) 7 A (l'2' ; 43) 



33'44' 



+ 



7 A (1'3'; 14) 7 A (2'4'; 23) - (1' <-> 2') - (1 <-> 2) + (1' <-> 2', 1 <-> 2)] } 



(10) 



and that 7 A=0 ° is equal to the bare anti-symmetrized two-particle interaction. In 
Eq. (10), we have already discarded a contribution containing the three-particle vertex 
7 A . This implements the above-mentioned truncation procedure leading to a closed set 
of differential equations for the self-energy and the two-particle vertex. It is strictly 
justified if the flowing couplings do not become too large (i.e., in general at small XT), 
since the function 7 A is zero at A = 00 and generated only by terms of third order in 
7 A . The truncation renders the FRG scheme an approximate method (containing at 
least all terms of order U 2 ), and one can obtain an approximation to the self-energy by 
integrating (in general numerically) Eqs. (7) and (10) from A = 00 to A = 0. 



3.2. Application to the SI AM 

In the following, we explicitly derive flow equations for the single impurity Anderson 
model defined in Sec. 2. To this end, it is useful to parametrize the two-particle vertex 
in a way which explicitly accounts for energy- and spin conservation 

7 A (a;^cri, u)' 2 o' 2 ; U\0-\,ix>iO-T) = 5t(u[ + oo' 2 — 0J1 — oo 2 ) 

(11) 

where St is the 5-function at zero temperature and the Kronecker symbol at finite T > 0, 
respectively. The quantities £/ A and £7^ denote three (new) independent functions. They 
fulfill certain symmetry relations which are discussed in detail in Appendix A. It will 
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prove useful both numerically and concept ionally (see Sees. 3.4 and 3.5 for details) to 
rewrite the arguments of U A in terms of three bosonic Matsubara frequencies (we skip 
the imaginary % for reasons of shortness) 

V\ = uj[ + uj' 2 , u 2 = w[ — uji, z/ 3 = uj' 2 — tO\. (12) 

The flow of the self-energy then takes the form [which follows from Eq. (7)] 

d A Z A (iu) = -Tj2e iQri [S A (iQ)U^{Q + uj,0,Q-uj) + S^{tQ)U^{Q + uj,0,±QTuj)} , 
in 

(13) 

where the upper sign holds for a =|, and we use the notation j =|. The initial condition 
reads T, A=OQ (iuj) = —U/2. The numerical integration, however, starts at some large but 
finite A . Due to the slow decay of the right-hand side of Eq. (13), the integration from 
A = oo to Ao does not vanish in the limit Ao — > oo, but rather tends to a finite constant. 
This constant can be computed analytically [16] and leads to the new initial condition 
E^°^°°(2u;) = 0. The flow equations for U A are determined by Eq. (10). The resulting 
expressions are lengthy and given explicitly in Appendix A. Together with Eq. (13), 
they form a closed set of differential equations, which is, however, still infinitely large. 
For any numerical treatment it is thus necessary to discretize the Matsubara frequency 
space using a finite number of points. This will be discussed in detail in Sec. 3.4. 




3. 3. Specification of a cutoff 

In a next step, it is necessary to specify the cutoff 0a- Since we are interested in treating 
finite temperatures T > 0, we choose a continuous function 

-ttT >\uj\- A 

-A) = <Jl/2-M^ -txT < \u\ - A < txT . (14) 

|w| - A > ttT 

Using this cutoff, the single-scale propagator Eq. (8) takes the form 

5»Jw -nT<\u\-A<nT ^ (lg) 
I otherwise 

and the Matsubara frequency sums appearing in Eqs. (13), (A.l), and (A. 3) reduce to 
two terms with frequencies uo for which — irT < \u\ — A < irT. 

In the zero-temperature limit T — >• 0, 0a becomes the usual step function (which 
is the cutoff that was previously employed for calculations carried out at T = [12, 16]), 
and one has to deal with products of ^-distributions S(\u \ — A) and functions / involving 
0(|cj| — A). These at first sight ambiguous expressions are well-defined and can be 
computed using Morris' Lemma [26], 

5 e (x - A)f [G e (x - A)] - S(x - A) f f(t) dt, (16) 
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where e is an arbitrary broadening parameter (in our case the temperature T), and 
5 f = Q' One obtains 



, A A t=o 5(\u\-A) 



S A (iu) = ^ '- , (17) 



and, defining 6(0) = 1/2, 



A A t=o 5(ki|-A) 0(|u; 2 |-A) 

S (lUJl)G {lU>2) 



[G°(iwi)] — E A (io;i) [G Q {iw2)]~ — £ A (iu;2) 

The ^-functions appearing in Eqs. (17) and (18) cancel the frequency integrals 
(originating from — > ^ f in the zero-temperature limit) on the right-hand side 
of the flow equations (13), (A.l), and (A. 3). 

The cutoff defined by Eq. (14) was successfully employed to tackle low-dimensional 
electron systems using the functional RG [12, 13, 14, 15, 16, 17] and is thus a very 
reasonable choice for the problem at hand. However, it is not always obvious that 
physical properties turn out to be independent of the actual realization of the cutoff 
procedure. Addressing this issue is subject of ongoing research. 



3.4. Numerical implementations 

Up to now, the only approximation involved was to discard the contribution of the 
three-particle vertex to the flow equation for 7 A , leading to the closed set of differential 
equations (13), (A.l), and (A. 3). This set is, however, infinitely large, so that for any 
numerical approach both S A and U A need to be parametrized by a finite number of 
flowing couplings. This introduces a second (numerical) approximation, and one has to 
carefully rule out the possibility of numerical artifacts. The present (technical) Section 
is devoted to this issue. 



3.4-1- Discretization of the frequency axis In this paper, we are mostly interested in 
the low-energy physics of the SIAM. Hence, it is reasonable to choose a parametrization 
of the Matsubara axis such that the low-frequency regime is better resolved than that of 
large \u\ (particularly having in mind the Kondo resonance at the chemical potential). 
At zero temperature, this is achieved by introducing a geometric mesh 

a n — 1 

lu = ±lu - — -, n=l...N, (19) 

with free parameters Uq > 0, a > 1, and N. At T > 0, the Matsubara frequency 
space is discrete by itself. The low-energy regime is then resolved by first accounting 
for the smallest iVo frequencies {dj ,^i, • • • ,^n }, next for max{iVo — S, 1} frequencies 
{ujno+a, &n +2A, ■ ■ ■ 1 uj n +(n -s)a} where A — 1 are left out in between each pair of kept 
frequencies, until the total number of frequencies if equal to N. Thus, at T > the 
imaginary axis is parametrized by four integer numbers N, N , S, and A. In both 
cases (T = and T > 0), on the right-hand side of the flow equations (13), (A.l), and 
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Figure 1. (Color online) (a) Imaginary part of the SIAM self-energy at U/T = 4 
and T = B = e = obtained from the FRG flow equations (13), (A.l), and (A. 3) for 
different discretization parameters N, w m i n , and u max . The two-particle vertex was 
parametrized using either the bosonic frequencies v [Eq. (12)] or the original arguments 
ui[, u' 2) and u>\ [Eq. (11)]. One should note that for the latter N = 50 roughly 
corresponds (concerning the computational effort) to N = 100 bosonic frequencies 
due to the lack of symmetries. Away from the discretization points determined by 
Eq. (19), both E A and 

J2 were evaluated using linear interpolation (LI) or the best 
matching frequency (BF). Note that the Bethe ansatz Kondo scale Eq. (27) is given 
by Tk/T = 0.29. FRG results were obtained using the modification introduced in 
Sec. 3.4.3. (b) The same, but comparing the effective mass m*(U) defined by Eq. (28) 
for different discretization parameters. For reference, the (numerically exact) result 
extracted from the numerical renormalization group is shown as well (thin dashed 
line). 



(A. 3) the two-particle vertex needs to be evaluated for arbitrary arguments which in 
general do not coincide with a particular point of the discrete frequency mesh. Hence, 
an interpolation procedure has to be devised. The set of flow equations then closes and 
can be solved numerically using standard Runge-Kutta routines. Simple interpolation 
procedures are to evaluate S A and U A at the nearest frequency or to interpolate linearly. 
For not too large U (roughly U/T < 7 for the problem at hand) one can verify that 
physical properties are equal for both of these choices, provided that the number of 
frequencies N is large enough [see Fig. 1(a)]. However, one observes that N has to be 
chosen much larger if no interpolation is performed, rendering it impossible to obtain 
convergent results (with respect to N, ujq and a) for larger U /T in that case. It is 
nevertheless reasonable to assume that this is just an issue of numerical resources. We 
thus refrain from implementing more complicated (e.g., spline) interpolation routines, 
and the data presented in this paper is obtained using linear interpolation of both S A 
and 7 A . 

As for the choice of the interpolation procedure it is of course imperative to check 
that physical results are not affected by the actual realization of the frequency mesh. 
In order to do so, we choose the parameters a>o and a (illustrating the procedure for 
T = 0) such that the smallest (largest) frequency cj min (u; max ) is at least 2-3 orders 
of magnitude smaller (larger) than the smallest (largest) energy scale of the system 
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(including the Kondo scale). We then increase the total number of frequencies N until 
physical properties remain unchanged up to a certain accuracy (usually one percent). 
Finally, we check that the same holds if u> m i n (c<j max ) is decreased (increased) further 
and if another kind of frequency mesh (e.g, a logarithmic mesh) is employed. As an 
example, the low-energy behavior of the FRG self-energy for different realizations of the 
Matsubara frequency discretization is shown in Fig. 1. 

As a final point, we ruled out the possibility that information about the system 
may be lost during the integration (and that this does not manifest at A = 0), or 
that the very choice of some exponentially-like frequency mesh is inappropriate even 
if one aims only at the low-energy physics. To this end, we implemented a frequency 
mesh which dynamically chooses its discretization points at each Runge-Kutta step 
and which adjusts the overall number of frequencies in a stepsize control manner such 
that a desired accuracy (e.g., for the self-energy) is constantly maintained during the 
integration process. We exemplary checked that our results are not changed by using 
these more elaborate (but slower) numerics. 

3.4-2. Parametrization of the two-particle vertex By discretizing the frequency axis, 
the two-particle vertex is parametrized by a three-dimensional set of points. However, 
the function 72 can be expressed in terms of three arbitrary independent arguments. 
The bosonic frequencies Eq. (12) are a very natural choice (see Sec. 3.5 for details). In 
addition, they are numerically favorable as all symmetries of 72 (see Appendix A) are 
preserved automatically (being sign change instead of interchange symmetries), which 
reduces the number of independent flow equations and the effort in finding nearest 
discretization points significantly. However, regarding 72 as a function of three of its 
original frequencies to is also an obvious possibility. Despite the lack of any fundamental 
arguments and the fact that symmetries are not preserved numerically, we exemplary 
investigated whether physical properties [particularly the low-energy behavior of the 
self-energy] are unchanged if 72 is parametrized in this way. We observe that for 
small to intermediate U/Y < 5, both parametrizations indeed give coinciding results. 
However, the number of discretization points N has to be significantly larger if the 
original frequencies uj instead of the bosonic frequencies are employed [compare the 
poor agreement between the solid and the long-dashed line in the main part of Fig. 1(a) 
with the perfect agreement even for \u\ > 0.02r between the two curves shown in 
the inset], supporting the assessment that the latter are the natural arguments of the 
two-particle vertex. In addition, it proves practically impossible to obtain coinciding 
results at larger U /Y. Even though it is again reasonable to assume that this is merely 
a question of numerical resources, one cannot finally verify that the strong-coupling 
results published in this paper are indeed independent of the actual parametrization of 
7^ (i.e., free of numerical artifacts). However, there are strong arguments favoring the 
use of the bosonic frequencies v. Thus, we stick to this choice, carefully ensuring that 
physical properties are independent of the corresponding parameters N, kj min and uj max 
[see Fig. 1(b) for the strong-coupling behavior of the effective mass]. 
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3.4-3. Modification of the flow equation for U A If one pursues the course of action 
outlined in Sec. 3.4.1, one can reliably extract physical properties of the SIAM (e.g., 
spectral functions) provided that the two-particle interaction is not too large (see Sec. 4). 
However, it turns out that if U/Y > 5, numerical integration of the flow equations fails. 
One observes that as A approaches some Ai ([/"), the size of integration steps needs 
to be continuously decreased in order to maintain a desired accuracy, until finally the 
machine precision is reached and no further progress can be made. This holds for 
various implementations of Runge-Kutta advancer routines [27]. However, one has to 
keep in mind that numerically one does not address the exact (up to truncation) flow 
equations (13), (A.l), and (A. 3), but an approximation induced by discretization of 
the frequency axis. Thus, before making any final statement one has to rule out the 
possibility that the breakdown of numerics is merely an artifact of this approximation. 
In order to do so, we repeated each Runge Kutta step with an increased number of total 
frequencies N until a desired accuracy for the self-energy was reached (this procedure 
was already mentioned in Sec. 3.4.1). It turned out that as A approaches Ai, N needs 
to be increased continuously until finally numerical resources are exhausted. Hence, it is 
impossible to say whether integration of the flow equations fails because of fundamental 
reasons (i.e., truncation of the infinite FRG hierarchy after second order) or because of 
the discretization of the frequency axis. 

Recently, a modification of the truncated FRG flow equations motivated by 
considering the fulfillment of Ward identities was suggested [28]. If (as in our case) 
the flow of the three-particle vertex is discarded, it consists of replacing 

S A _> _dG^ = s a _ g^gA (20) 
dA dA 

on the right-hand side of the flow equation for the two-particle vertex. Obviously, the 
additional term proportional to Q a Tj A Q a is at least of third order in the bare interaction 
strength. Indeed, implementing the flow equations (A.l) and (A. 3) for U A with S A 
replaced by —d\Q A gives results which coincide [on the scale of the corresponding Figure 
such as Figs. 2(a) or 3(b)] with the original ones for U/Y < 3 [29]. One observes, however, 
that numerical integration of the modified flow equations is possible even if U/Y > 5 
and gives results in better agreement with the NRG reference for U/Y ~ 3 — 5. Hence, 
we pragmatically employ these equations for the rest of the paper, keeping in mind that 
it is yet unclear why the original scheme breaks down for U/Y > 5. 

3.5. Further approximations 

As mentioned in the previous Section, the number of total discretization points N has 
to be chosen such that physical properties remain unchanged if N is increased further. 
For the problem at hand, it turns out that this holds if N « 50. Hence, one has to 
numerically integrate approximately 10 5 coupled ordinary differential equations, each 
containing an integral on the right-hand side. The latter originates from the second 
term in Eq. (20). Using parallelized code running on typically 8 CPU cores, numerics 
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needs a couple of days of computer time for a single set of parameters. Thus, it is 
desirable to devise further approximations which allow a fast qualitative overview of the 
general physics. One such approximation (which is to be called appr. 1) can be obtained 
by setting 

Viti = (21) 

in the terms proportional to V A (iuJ, ± ioj) on the right-hand side of Eqs. (A.l) and 
(A. 3). The two-particle vertex U A is then no longer parametrized by three independent 
arguments but by three functions each depending on one variable This reduces 
the numerical effort dramatically. A motivation for this approximation can be given 
by noting that the three Feynman diagrams [corresponding to the particle-particle and 
the two particle-hole - terms on the right-hand side of Eqs. (A.l) and (A. 3)] which 
contribute to 72 in second-order perturbation theory each depend on a single frequency 
^'=1,2,3 only. Appr. 1 is thus justified for small U. Indeed, it turns out that for U/T < 3, 
one can obtain quantitatively reliable results within minutes of single-core CPU time (see 
Sec. 4.1.3). In addition, the connection to perturbation theory a posteriori illustrates 
that the bosonic frequencies z/j [Eq. (12)] can be regarded as the natural arguments of 
the two-particle vertex 7^. 

An even simpler set of flow equations (which is to be called appr. 2) can be obtained 
by completely disregarding the frequency- dependence of 7^ (up to energy conservation). 
For the problem at hand, evaluating Eqs. (13), (A.l), and (A. 3) for zero frequency yields 
(see also Refs. [16, 30]) 

d A Z A = -TU^^) , 

iul (22) 

iuj a 

Thus, the flowing self-energy remains frequency-independent, and S A=0 can be regarded 
as an effective-single particle potential. The spectral function is then by construction 
a Lorentzian of width 2r for arbitrary U, illustrating that one cannot reliably extract 
finite-frequency properties from a frequency-independent FRG truncation scheme. On 
the other hand, Kondo pinning of the spectral weight at the chemical potential is 
described well by this approximation (even though this is a strong coupling effect), and 
the linear-response conductance as a function of the impurity energy e was computed 
accurately for a variety of quantum dot systems [15, 16, 17]. 

4. Results 

In this Section, we present results for physical properties of the SIAM which we 
obtain from numerically integrating the flow equations (13), (A.l), and (A. 3) with 
the modification introduced in Sec. 3.4.3. Our findings are benchmarked against data 
obtained from the framework of the numerical renormalization group. First, we discuss 
the spectral function (Sec. 4.1) as well as the average impurity occupation number and 
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Figure 2. (Color online) (a) Impurity spectral function p(u>) of the SIAM for weak 
Coulomb repulsion U/T = 2, zero temperature and zero magnetic field at the particle- 
hole symmetric point and for finite impurity energy e = T. FRG results were obtained 
using Pade approximation of the Green function calculated from Eqs. (13), (A.l), 
and (A. 3). The discretization parameters read N = 75, u) m i n /T = 1 • 10~ 5 , and 
^max/r = 2259. For reference, spectral functions calculated from the restricted 
Hartree-Fock approximation (RHF) are shown as well [33]. (b) The same, but for 
intermediate coupling U/T = 4 and both zero and finite temperature T. The frequency 
discretization for T > is determined by N = 81, Nq = 18, S = 2, and A = 2. 

the zero-temperature spectral weight at the chemical potential (Sec. 4.2), in particular 
in the regime of small to intermediate U/T < 5. Then, we investigate whether the 
FRG approximation contains an exponential energy scale in the strong coupling limit 
(Sec. 4.3). In both cases, we elaborate on how our results are modified by the additional 
approximations introduced in the previous Section. 

4-1. Spectral functions 

4-1.1. Pade approximation By construction, our functional RG scheme gives the self- 
energy Y>(iu) and the interacting Green function Q{iu) in Matsubara frequency space. 
Hence, an analytic continuation needs to be performed in order to obtain the spectral 
function 

p(u) = — Jmg(u + irj). (23) 

7T 

Since Q{iuj) is not known analytically but only for a certain discrete number of points, 
such a continuation has to be carried out numerically. In this paper, we compute p(u>) 
using Pade approximation [31, 32]. We address the possibility of numerical artifacts 
(analytic continuation being an ill-posed problem) by computing each particular curve 
for various (at least two) realizations of the Matsubara frequency discretization. In 
addition, we check that calculating p from Pade approximation of the Green function 
Q{iuj) and the self-energy E(iu;) (using an odd and even number of discretization points, 
respectively) gives coinciding results. In short, analytic continuation turns out to be 
stable at zero temperature, small U/T < 5, small e < U/2, and small B < U/10, whereas 
it is more difficult to reliably calculate spectral functions for larger Coulomb interaction 




Figure 3. (Color online) The same as Fig. 2, but additionally showing spectral 
functions obtained from the FRG approximation 1 introduced in Sec. 3.5. 



and particularly finite T > 0. 

4-1-2. Spectral functions In Figs. 2, 3, and 4, we compare p(u) obtained from the NRG 
framework with FRG results. For reference, data calculated using the restricted Hartree- 
Fock approximation [33] is shown as well. For small Coulomb interaction U/Y = 2—4 the 
NRG and frequency- dependent FRG schemes agree perfectly. This holds particularly if 
both the impurity energy and the magnetic field are not too large [see Fig. 2(a) and 3 
for finite e and Fig. 4(b) for finite E\. For U/Y > 3, the spectral functions from both 
schemes begin to deviate if either e/Y or B/Y is increased sizably [see Fig. 3(b) for 
the e-dependence] . However, at the same time the Pade approximation is observed to 
become unstable, rendering it reasonable to consider quantities which can be calculated 
directly from the imaginary axis (such as the zero-temperature spectral weight at the 
chemical potential and the average occupation number) in order to assess the accuracy 
of the FRG approximation for large e/Y and B/Y. This is done in Sec. 4.2. 

It is particularly important to point out that the FRG approximation yields 
accurate results both for T = and finite temperatures [see Figs. 2(b) and 3(a) for 
T > 0] . One should note that accounting for the frequency-dependence of 7^ within the 
FRG truncation scheme is imperative in order to properly extract a spectral function at 
T > (the latter being a Lorentzian of constant height 1/tt in the frequency-independent 
appr. 2 of Sec. 3.5). In that respect, it would be desirable to further investigate 
the temperature-dependence of physical properties (e.g., Fermi-liquid behavior). The 
Pade approximation, however, turns out to be particularly ill-controlled for T > 
and extracting p(uS) for a variety of parameters is practically impossible. In contrast, 
computing the NRG Matsubara Green function from the real-axis data (by virtue of 
a Hilbert transformation) is numerically well-controlled, and at small to intermediate 
U/Y < 4 perfect agreement with FRG results is observed for arbitrary temperatures. 

For some larger U/Y = 5, FRG data is still in good agreement with the NRG 
reference [see Fig. 4(a)]. In contrast, the spectral function obtained from second- 
order perturbation theory [34] already deviates sizably from the NRG curve, and one 
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Figure 4. (Color online) (a) The same as Fig. 3, but for larger U/T = 5. The second- 
order perturbation theory result is shown as well (PT). (b) Spectral functions obtained 
from appr. 1 in presence of a magnetic field B. The upper panel shows the spin-up 
component P](lo) only. 



observes that this disagreement becomes even worse if e is shifted away from particle- 
hole symmetry. If U is increased further towards the strong-coupling limit, we observe 
that even &te = B = T = the Pade approximation becomes unstable for large energies 
\u\ > 2T K (T K being the Kondo scale). Hence, it is impossible to address the question 
whether the FRG scheme describes Hubbard bands, whereas a resonance at the chemical 
potential can still be extracted reliably. As mentioned above, the agreement between the 
FRG and NRG frameworks deteriorates for U/Y > 5. This will be illustrated explicitly 
(also in relation to perturbation theory) in the next Section where the width of the 
oj = (Kondo) - resonance is compared. In short, one can state that for the problem at 
hand the FRG approximation is accurate for small to intermediate Coulomb interaction 
U/T < 5. 

As mentioned in Sec. 3.4.3, performing FRG calculations without replacing the 
single-scale propagator by Eq. (20) gives results which almost coincide with those 
presented in Figs. 2(a) and 3(b). The deviation between both schemes is larger (though 
still not sizable) in Fig. 2(b) (where U/Y — 4), and the data computed with S A — > d\Q A 
turns out to be quantitatively better. 

4-1.3. Spectral functions - appr. 1 As mentioned above, producing one of the curves 
shown in Fig. 2 usually takes a time span which is of the order of days. Hence, 
it is desirable to determine the accuracy of the (numerically far less demanding) 
approximations devised in Sec. 3.5. This is done in Figs. 3 and 4. One observes that even 
for U/Y = 4 and U/Y = 5, appr. 1 gives results which agree well with those obtained 
from Eqs. (13), (A.l), and (A. 3) and with NRG, respectively. This holds for arbitrary 
impurity energies e and both zero- and finite temperatures T. If U/Y < 3, the spectral 
functions computed from the NRG, FRG and FRG appr. 1 become indistinguishable. 
Hence, the latter provides an efficient tool to extract p(oj) provided that the Coulomb 
interaction is not too large (U/Y < 5 in case of the SIAM). 
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In contrast, the spectral function is by construction a Lorentzian of interaction- 
independent width 2T and temperature-independent height 1/(ttT) if one completely 
discards the frequency-dependence of the two-particle vertex and the self-energy (appr. 2 
introduced Sec. 3.5). Hence, this approximation is not suited to compute finite-frequency 
properties of the SIAM. On the other hand, it was observed in a previous publication 
(Ref. [16]) that for T = the spectral weight at to = is described well by the 
frequency-independent FRG scheme even in the strong-coupling limit, and the zero- 
temperature linear-response conductance G(e) (which is a zero-energy property) shows 
the characteristic Kondo plateau in quantitative agreement with NRG calculations [see 
the upper panel of Fig. 5(b) for some fairly small (in the context of strong-coupling 
physics) U/T — 6]. 

4-2. Spectral weight at uj = and average occupation 

In Fig. 5, we show the zero-temperature spectral weight at the chemical potential 
p(u = 0) and the average occupation number (n) as a function of both the magnetic 
field and the impurity energy. We exclusively employ the numerically less demanding 
FRG appr. 1 in order to tackle the large number of parameters but exemplary checked 
that the full (frequency-dependent) approximation gives coinciding results (in agreement 
with the statements of Sec. 4.1.3). Both p(0) and (n) can be extracted directly from the 
imaginary axis and can thus be computed without the need for an ill-controlled analytic 
continuation. In particular, the average occupation can by definition be obtained by 
integrating the Matsubara Green function over the imaginary axis (focusing exclusively 
on T = 0): 

M = 7T f e^Mdcu . (24) 



2tt _ 

In comparison with NRG as well as with exact Bethe ansatz data (taken from Ref. [36]), 
one observes that the frequency-dependent FRG scheme describes both p(u = 0) and 
(n) very accurately for small U/T = 2 and intermediate U/T = 4 — 6. 

It is also possible to extract the average occupation (n) from a Friedel sum rule [5], 



11 

("Wfsr = o arctan 



e + aB/2 + Re E(iri) 



(25) 



as well as from the derivative of the grand canonical potential Q with respect to the 
single-particle energy e [37]: 

(n,+n l ) n = ^ ■ (26) 

Since the truncated FRG is not a conserving approximation, these different ways to 
compute the average occupation do not necessarily give coinciding results. However, 
using FRG appr. 1 one observes that both (n) and (n)n agree well with the Bethe 
ansatz result even for U/T = 6 [see Fig. 5(b)]. In contrast, (n)n differs sizably from the 
exact data if the simpler frequency-independent FRG scheme (appr. 2) is employed. 
This illustrates that the FRG approximation to the grand canonical potential is 
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Figure 5. (Color online) (a) Comparison of NRG and FRG data for the zero- 
temperature spectral weight p{u> = 0) and for the average occupation number (n) 
[determined by Eq. (24)] as a function of the magnetic field B for different parameters, 
(b) Average occupation extracted from integrating the Matsubara Green function 
[Eq. (24)], from a Friedel sum rule [Eq. (25)], and from the derivative of the free 
energy [Eq. (26)] in comparison with Bethe ansatz data taken from Ref. [36]. FRG 
calculations were carried out using the approximation 1 as well as the frequency- 
independent scheme (appr. 2). 



systematically improved by accounting for the frequency dependence of the two-particle 
vertex. In contrast, the Friedel sum rule (n) = (?t.)fsr is fulfilled analytically within the 
(effectively noninteracting) approximation 2, whereas the frequency-dependent schemes 
yield a self-energy Yi(irj) whose real part deviates from the Bethe ansatz result (at 
fairly large U/T = 6) for intermediate e/T. This agrees with the observation that (at 
U/T > 3) the FRG spectral function is more accurate close to particle-hole symmetry 
than for larger e. 



4-3. The Kondo scale 

In this Section, we investigate whether the FRG framework contains an exponentially 
small energy scale governing the low-energy physics in the strong-coupling limit. From 
the exact Bethe ansatz solution of the SIAM it is known that this so-called Kondo 
temperature T K is given by [5] 

T K = v^7lV2exp \U 2 -4e 2 |] . (27) 

In Sec. 4.1 we have already noted that NRG and FRG results begin to deviate seriously 
for U/T > 5. However, it was observed in a previous publication (Ref. [16]) that an 
exponential energy scale can be extracted out of zero-energy properties (e.g., the static 
spin susceptibility) using a frequency- independent FRG truncation scheme (appr. 2). 
Hence, it is reasonable to address the same issue using the frequency-dependent FRG 
approach introduced in this paper. 
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Figure 6. (Color online) (a) Comparison of NRG and FRG results (at e = T = B = 0) 
for the effective mass m*(U) defined by Eq. (28). In the strong coupling limit 
[//r> 1, the NRG scheme gives m* oc , whereas the FRG approach does not show 
exponential behavior. The FRG results were obtained using discretization parameters 
N = 75, cj m j n /r = I • 10~ 5 , and uJ max /T = 2259 [38]. The effective mass extracted 
from second-order perturbation theory is shown as well (PT). (b) The same, but for 
the static spin susceptibility defined by Eq. (29). FRG results were computed using 
N = 50, ui min /T = 1 • I0~ 5 , and w max /r = 980. Data obtained from the frequency- 
independent truncation scheme (appr. 2) is shown as well. For reasons of convenience, 
symbols (connecting lines) are not shown in the main part (the inset). 



4-3.1. Effective mass A characteristic low-energy property governed by the Kondo 
temperature is the width of the zero-frequency resonance in the spectral function at 
e = T = B = 0. Equivalently, one can consider the effective mass m*(U) defined by 

d[lmE(iu)] d[lmg(iu)} 
m {U) = 1 .=o + = fa - • {2t 



duj 



uj=0+ 



This quantity can be accessed directly within the Matsubara FRG formalism, whereas 
it is computed from the NRG spectral functions using a Hilbert transformation. The 
latter is much more controlled than analytic continuation from the imaginary to the 
real axis, so that numerical artifacts originating from Pade approximation are ruled out 
from the beginning. The results are shown in Fig. 6(a). One observes that m*(U) shows 
quadratic behavior for small U, and that data obtained from the NRG, FRG and FRG 
appr. 1 agree well for U/T < 5. In contrast, there is already a sizable deviation between 
the NRG reference and perturbation theory [34] for U/T >3. In the Kondo regime, the 
effective mass computed from the NRG approach is inversely proportional to the Kondo 
temperature. In contrast, the FRG scheme does not produce an exponential behavior. 
It deviates seriously from the NRG reference for U/T > 5, a mere observation being 
that the appr. 1 gives results which are quantitatively better. The same statements 
hold if one considers the width of the spectral function itself. It turns out that in 
the strong coupling limit the latter can be obtained by scaling l/m*(U) with a U- 
independent factor. Recalling the observation that in general the Pade approximation 
becomes unstable for large U, this is a clear indication that it is yet possible to reliably 
extract a zero-frequency resonance from Matsubara FRG calculations. 
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4-3.2. Static spin susceptibility A second quantity governed by the Kondo energy scale 
is the static spin susceptibility x- It is defined as 

x(.u) = -liM-iM> 

where (n a ) denotes the average occupation number of electrons with spin direction a at 
the impurity [extracted using Eq. (24)]. The results are shown in Fig. 6(b). Again, one 
observes good agreement between NRG and FRG data for U/F < 5. The susceptibility 
derived from the NRG approach becomes proportional to the inverse Kondo temperature 
in the strong-coupling limit, whereas the FRG framework does not yield exponential 
behavior. In contrast, the susceptibility computed using the frequency-independent 
FRG scheme (appr. 2) is of the form exp(— cU/F) with c~l/7r^7r/8. The Kondo scale 
is thus contained at this simple level of truncation but not within the (by construction 
'higher-order') frequency- dependent approximation. On the other hand, one observes 
that at small U/F < 5, the agreement with the NRG reference improves significantly 
by employing the latter [see the inset of Fig. 6(b)]. The same holds if within appr. 2 
the Kondo temperature is defined as the magnetic field Bk necessary to suppress the 
spectral weight at the chemical potential down to l/(27rr). Again, it turns out that 
Bk oc exp[— U / (7rr)] (see Ref. [16]), but this exponential behavior is lost if the frequency- 
dependence of 7^ is accounted for. 

5. Conclusions and Outlook 

In this paper, we have introduced a truncation scheme for the infinite hierarchy of FRG 
flow equations which accounts for the frequency dependence of the two-particle vertex. 
The flow equations for the single impurity Anderson model were derived explicitly. Using 
a variety of parametrization procedures, we carefully addressed the issue of numerical 
artifacts originating from the discretization of the Matsubara axis. We showed that at 
intermediate Coulomb interaction U/F < 5 one can obtain data which is independent of 
any numerical parameters. In contrast, different ways to parametrize the three frequency 
arguments of the two-particle vertex do not give coinciding results in the strong-coupling 
regime, even though it is reasonable to assume that this is merely a question of numerical 
resources. In addition, there are strong (conceptional and practical) arguments favoring 
the use of a certain set of bosonic frequencies. We thus employed this parametrization 
to carry out calculations at large U, carefully ensuring that the results are independent 
of any remaining numerical parameters (particularly the number of frequencies N). 

In general, it turned out that at small to intermediate U/F < 5 the FRG 
approximation works well (benchmarking our results against numerical RG reference 
data), and properties such as spectral functions, which are certainly of both 
experimental and theoretical interest, can be computed accurately for arbitrary 
parameters (particularly finite temperatures). Using a simplified set of flow equations, 
such reliable calculations are possible within minutes of CPU time. In contrast, 



(29) 
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it proved impossible to tackle the strong-coupling limit, and zero-energy aspects of 
Kondo physics which are captured by a simple frequency-independent FRG scheme 
(e.g., the static spin susceptibility x being governed by an exponential energy scale 
[16]) are no longer described by the 'higher-order' frequency-dependent approximation. 
However, it is imperative to keep in mind that it was neither possible to clarify why 
one particular (unmodified) FRG scheme breaks down in the strong-coupling limit 
nor practically manageable to obtain data independent of all numerical parameters 
(namely the parametrization of the two-particle vertex). On the other hand, at small to 
intermediate U the agreement with NRG or Bethe ansatz data is improved quantitatively 
by employing the more elaborate finite-frequency scheme. These observations are 
consistent with the fact that the latter contains all terms up to order U 2 (but 
gives results superior to second-order perturbation theory), whereas the frequency- 
independent approximation is only correct to first order in U. However, the frequency- 
independent approach does not suffer from typical mean- field artifacts (e.g., breaking of 
spin symmetry) and gives results which are in significantly better agreement (compared 
to Hartree-Fock) with reference calculations (concerning the low-energy physics of 
quantum dots [39]). It can thus be pragmatically viewed as a reliable tool to derive 
effective noninteracting parameters which accurately describe zero-energy aspects (e.g., 
the T = linear-response conductance) of correlation phenomena (RG enhanced 
Hartree-Fock theory). 

Concerning prospects for future work, it would be interesting to apply the 
frequency-dependent FRG scheme presented in this paper to such multi-impurity 
systems which cannot be accessed easily using the NRG framework (the numerical effort 
growing only as a power law but not exponentially with the number of impurities and 
the number of channels). In the context of quantum wires, a renormalization group- 
based approach to any microscopic model is inherently necessary because of infrared 
divergences in perturbation theory. Generalizing the method to the non-equilibrium 
situation is certainly an issue of interest. There is ongoing work in all these directions 
(see Refs. [40, 41, 42]). From a conceptional point of view, it would be favorable to 
set up an FRG scheme in real frequency space which does not suffer from the need for 
an ill-controlled (particularly at finite temperatures) analytic continuation. There is 
ongoing research in this direction as well, and first results were published [40, 42]. 
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Appendix A. Flow equation for U A 

The flow of the three independent parts which parametrize the two-particle vertex 7 A of 
the single impurity Anderson model is determined by Eq. (10). For U^ = ^ j_ one obtains 

= - T ^2[ 'P£ a {iui,ivi-iu>)U£{vi,ui + v_ ++ ,v +++ -u>) U£(v lt u) + !/_+_, ui + v__ + ) 
+ { Ve (T {iu),iv2+iu)U£{u) + v ++ _,i>2,u + V- ++ ) U^(ui + v +++ , -u 2 , ui + v_ + _) 
+V^(iu), iv 2 + iui) Ufl(u + ±i/ 2 ,±w ± U^(u> + v +++ , T" 2 , ±ui ± 

-(l/ 2 -> -I> 2 , V3 -* -^3) - [(^2 *-> ^3) ~ (^2 -> ~V 2l U 3 -> -I/3)] } 

(A.l) 

and the initial condition reads /7 A=00 (z/i, z/2, z/3) = U A °^°°{vi, z/ 2 , z/3) = 0. Again, the 
upper sign holds for the spin up component. We have introduced shorthands like 
z/-+- = [—v\ + z/ 2 — Vs)/2 and defined the quantity 

V^iu^) = S^iujg^iu*). (A.2) 
The flow equation for has the form 
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(A.3) 

with the initial condition being t7^ =<X) (z/i, z/ 2 , z/3) = i7-j^ 0_>oo (z/i, z/ 2 , z/3) = U. The functions 
f/ A obey the following symmetry relations (reflecting the symmetries of 7 A ): 

t/ CT Vl,^ 3 ) = ^(^,±^,±^3) = [U^-V!,!^,!*)]* = -UZ{v u v 3 ,V2) (A.4) 
for cr =|, I, and 

Ufax, u 2 , u 3 ) = Ufifa, -u 2 , u 3 ) = [Uf^-uu u 2 , -u 3 )]\ (A.5) 

These relations can be exploited in numerics and hold in the most general case. For 
many parameters of interest, however, additional symmetries are fulfilled. In particular, 
for particle-hole symmetry e = 0, 

[tff fa, z/ 2 , vs)] * = C/f z/ 2 , z/ 3 ), f/^z/x, z/ 2) z/ 3 ) = f/u(-z/i, z/ 2 , z/ 3 ), (A.6) 
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and the self-energy components are related by —T,*(iu) = Ej(io;). In absence of a 
magnetic field, spin symmetry leads to 

Uf(v u vi, i*0 = Uffa, »2, va) = Utlfa, V2, u 3 ) - U^{y x , u 3 , v 2 ), (A.7) 

and 

Uk(u 1 ,v 2 ,v 3 ) = Ufi(v 1 ,v 2 ,-v 3 ). (A.8) 

If in addition e = holds, the function is purely real (by virtue of the symmetry 
under v\ — > — V\), whereas the self-energy is purely imaginary. This fact can be 
conveniently exploited to speed up numerics. 
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